/*
TRABALHO FINAL DE PROGRAMACAO CONCORRENTE 2013
FEITO POR:
Danilo Augusto Silva, 6418678
Felipe Fukakusa, 6426976
Jader Bruno Pereira Lima, 6792346
*/

#include <stdlib.h>
#include <stdio.h> 
#include <string.h> 
#include <math.h>
#include <sys/time.h>

//=========FUNCOES===========//

//TestarConvergencia: Testa se a matriz A fornecida é viavel para que se convirja para
//um resultado pelo método de Jacobi-Richardson
bool TestarConvergencia(float*);
//Teste de verificação da linha passada no arquivo (como dado no enunciado)
void RowTest(float *resultado, float *A, int tamanho);

//========VARIAVES GLOBAIS=========//
//ORDEM DA MATRIZ A
int J_ORDER;
//LINHA A SER TESTADA
int J_ROW_TEST;
//NUMERO MAXIMO DE ITERACOES
int J_ITE_MAX;
//ERRO MAXIMO PERMITIDO
float J_ERROR;

//========FUNCOES EM CUDA=========//

//JacobiRichardsonGPU: Realiza a teracoes do método de Jacobi-Richardson para cada indice em cada nó da
//placa gráfica
__global__ void JacobiRichardsonGPU(float* A, float* x, float* b, int J_ORDER, float J_ERROR, int J_ITE_MAX){
    //=======VARIAVEIS=======//
    float sum = 0;  //Auxiliar para a soma no método
    float temp;     //Auxiliar para o método
    float error;    //Calculo do erro a cada iteracao
    int it = 0;

    //RECUPERA O INDICE A SER PROCESSADO
    int indice = threadIdx.x + blockDim.x * blockIdx.x;
    
    //LAÇO PRINCIPAL, FAZER ATÉ que a diferença entre os valores de iterações vizinhas seja menor que o erro
    do{
        it++;
        sum = 0;
        for (int k = 0; k < J_ORDER; k++){
            if(indice!=k)
                sum += A[(indice*J_ORDER)+k] * x[k];
        }
        
        temp = (b[indice] - sum)/A[(indice*J_ORDER) + indice];
        
        //calculo do erro
        error = fabs(x[indice] - temp);
        
        x[indice] = temp;
    } while(error > J_ERROR && it <= J_ITE_MAX);
    
}

//FUNCAO PRINCIPAL
int main(int argc, char** argv) {
    //========VARIAVEIS=========//
    int SIZE_A, SIZE_MEMORY_A;  //TAMANHO DAS MATRIZES DO SISTEMA
    int SIZE_X, SIZE_MEMORY_X;
    int SIZE_B, SIZE_MEMORY_B;
    int SIZE_Y, SIZE_MEMORY_Y;
    
    int i, j, k; //CONTADORAS
    
    struct timeval startTime;   //VARIAVEIS DE TEMPO
	struct timeval endTime;
    struct timeval startGPU;
    struct timeval endGPU;
    
    float *A, *B, *X, *Y_GPU;   //MATRIZES DO SISTEMA

    FILE *file;

    const int max_threads_por_bloco = 200;  //NUMERO DE THREADS POR BLOCO NA PLACA GRAFICA
    
    //INICIA CONTAGEM DO TEMPO TOTAL
    gettimeofday( &startTime, NULL );
    double start = startTime.tv_sec + ( startTime.tv_usec / 1000000.0 );
    
    if(file = fopen(argv[1], "r")){
                
        //LEITURA DOS DADOS DO ARQUIVO
        fscanf(file, "%d\n", &J_ORDER);
        fscanf(file, "%d\n", &J_ROW_TEST);
        fscanf(file, "%f\n", &J_ERROR);
        fscanf(file, "%d\n", &J_ITE_MAX);
        
        printf("J_ORDER = %d\n", J_ORDER);
        printf("J_ROW_TEST = %d\n", J_ROW_TEST);
        printf("J_ERROR = %f\n", J_ERROR);
        printf("J_ITE_MAX = %d\n", J_ITE_MAX);
        
        printf("\nEXECUTANDO PROGRAMA GAUSS_SEIDEL EM CUDA...\n");
        
        //CRIACAO DAS MATRIZES DO SISTEMA
        SIZE_A = J_ORDER * J_ORDER;
        SIZE_MEMORY_A = sizeof(float) * SIZE_A;
        A = (float*) malloc(SIZE_MEMORY_A);
    
        SIZE_X = J_ORDER;
        SIZE_MEMORY_X = sizeof(float) * SIZE_X;
        X = (float*) malloc(SIZE_MEMORY_X);
        
        SIZE_B = J_ORDER;
        SIZE_MEMORY_B = sizeof(float) * SIZE_B;
        B = (float*) malloc(SIZE_MEMORY_B);
    
        SIZE_Y = J_ORDER;
        SIZE_MEMORY_Y = sizeof(float) * SIZE_Y;
        Y_GPU = (float*) malloc(SIZE_MEMORY_Y);
        
        //LEITURA DA MATRIZ A
        k = 0;
        for(i = 0; i < J_ORDER; i++){
            for(j = 0; j < (J_ORDER - 1); j++){
                fscanf(file, "%f ", &A[k++]);
            }
            fscanf(file, "%f\n", &A[k++]);
        }
        //LEITURA DA MATRIZ B
        for(i = 0; i < SIZE_B; i++){
            fscanf(file, "%f\n", &B[i]);
        }
        //INICIALIZACAO DA MATRIZ X (INICIALAMENTE COM 0,0)
        for(i=0; i < SIZE_X;i++)
        {
            X[i]=0.0;
        }

   

        
        

        if(TestarConvergencia(A)) {
            //INICIALIZACAO DAS MATRIZES QUE SERAO USADAS NA PLACA
            float* A_GPU;
            cudaMalloc((void**) &A_GPU, SIZE_MEMORY_A);
        
            float* X_GPU;
            cudaMalloc((void**) &X_GPU, SIZE_MEMORY_X); 
    
            float* B_GPU;
            cudaMalloc((void**) &B_GPU, SIZE_MEMORY_B);
    
            
            //COPIA AS MATRIZES PARA AS MATRIZES DO DISPOSITIVO
            cudaMemcpy(A_GPU, A, SIZE_MEMORY_A, cudaMemcpyHostToDevice);
            cudaMemcpy(X_GPU, X, SIZE_MEMORY_X, cudaMemcpyHostToDevice);
            cudaMemcpy(B_GPU, B, SIZE_MEMORY_B, cudaMemcpyHostToDevice);
    
            //NUMERO DE BLOCOS DE PROCESSAMENTO NA PLACA
            int blocos = J_ORDER / max_threads_por_bloco;

            //INICIA CONTAGEM DO TEMPO DE PROCESSAMENTO NA PLACA
            gettimeofday( &startGPU, NULL );
            double startPlaca = startGPU.tv_sec + ( startGPU.tv_usec / 1000000.0 );
        
            //FUNCAO DE PROCESSAMENTO EM CUDA
            printf("\nEXECUTANDO PROCESSAMENTO NA PLACA GRAFICA...\n");
            JacobiRichardsonGPU<<< blocos, max_threads_por_bloco >>>(A_GPU, X_GPU, B_GPU, J_ORDER, J_ERROR, J_ITE_MAX);

            //FINALIZA A CONTAGEM DE TEMPO DE PROCESSAMENTO NA PLACA
            gettimeofday( &endGPU, NULL );
            double endPlaca = endGPU.tv_sec + ( endGPU.tv_usec / 1000000.0 );
    
            //COPIA AS MATRIZES DO DISPOSITIVO PARA O PROGRAMA
            cudaMemcpy(Y_GPU, X_GPU, SIZE_MEMORY_Y,cudaMemcpyDeviceToHost);
    
            //TESTA O RESULTADO DO SISTEMA
            RowTest(Y_GPU, A, J_ORDER);
    
            //FINALIZA A CONTAGEM DE TEMPO TOTAL
            gettimeofday( &endTime, NULL );
            double end = endTime.tv_sec + ( endTime.tv_usec / 1000000.0 );
    
        
            printf("B[ROW_TEST] = %f\n\n", B[J_ROW_TEST]);
        
            printf( "Tempo de execução total: %.6lf segundos\n", end - start );
            printf( "Tempo de execução de processamento de todos os nós na placa: %.6lf segundos\n", endPlaca - startPlaca );

            //DESALOCACAO DAS MATRIZES
            free(A); 
            free(B);
            free(X);
            cudaFree(A_GPU);
            cudaFree(B_GPU);
            cudaFree(X_GPU);
            cudaFree(Y_GPU);
        
            cudaThreadExit();

        }
        else{
            printf("Sistema nao converge");
        }
    }
    else{
        printf("ERRO ao abrir o arquivo\n");
    }
}

//TESTA SE A MATRIZ A CONVERGE PARA UMA SOLUCAO
bool TestarConvergencia(float *A){
    //VARIAVEIS AUXILIARES
    float linha_soma, coluna_soma;
    
    for(int i=0; i<J_ORDER; i++){
        linha_soma=0;
        
        for(int j=0; j<J_ORDER; j++)
            if(i != j)
                linha_soma += A[i*J_ORDER+j]; 
                
        linha_soma=linha_soma/A[i*J_ORDER+i];
        
        if(linha_soma>1) {
            for(i=0;i<J_ORDER;i++){ 
                coluna_soma=0;
                
                for(int j=0; j<J_ORDER; j++)
                    if(i != j)
                        coluna_soma += A[j*J_ORDER+i]; 
                
                coluna_soma=coluna_soma/A[i*J_ORDER+i];
                
                if(coluna_soma>1) 
                    return false;
            }
            return true;
        }
    }
    return true;
}

void RowTest(float *resultado, float *A, int tamanho){
    float result = 0;
    int i;
    
    for(i = 0; i < tamanho; i++){
        result += resultado[i]*A[J_ROW_TEST*tamanho + i];
    }
    printf("\nVERIFICACAO DO RESULTADO:\n");
    printf("ROW TEST: %f\n", result);
}

